!############################################################
 module flow
! Solution of the Stokes equation
!############################################################
   implicit none
   integer :: m1, ndim, nnz
   real (kind=8) :: loadx, loadz, etaN, etaS, etaE, etaW
   real (kind=8), dimension(:), allocatable :: rhs, sol, upom, vpom, ppom
   real (kind=8), dimension(:,:), allocatable :: epsII, epsII_previt, u_previt, v_previt, mtrx, b
   integer, dimension (:), allocatable :: rowind, columns, irn_m, jcn_m
   real (kind=8), dimension(:), allocatable :: nnzval
!###########################################################################

!*******************
  contains
!*******************

!============================================================
 subroutine stokes_alloc
! Allocate the variables needed to solve the Stokes equation 
!
! ndim => integer. Size of the vector of the unknowns
!                  nz*(nx+1) nodes for vertical flow
!                  nx*(nz+1) nodes for horizontal flow
!                  (nx+1)*(nz+1)-1 nodes for pressure 
! NB p is presicribed in the upper left corner to ensure uniqueness
!
! rhs:  right hand side
!============================================================ 
   use share; use domain
      
   ndim=3*nx*nz+2*(nx+nz)
   m1=3*nz+4 
   allocate(sol(ndim),rhs(ndim),b(ndim,1),alpha(0:nx+1,0:nz+1),kap(0:nx+1,0:nz+1),           &
            upom(nx*(nz+1)),vpom(nz*(nx+1)),ppom((nx+1)*(nz+1)-1),                           &
            u(0:nx+1,0:nz+1),v(0:nx+1,0:nz+1),p(nx+1,nz+1),urms(0:nx+1,0:nz+1),              &
            epsII(0:nx+1,0:nz+1),epsII_previt(0:nx+1,0:nz+1),visc_mech(0:nx+1,0:nz+1),       &
            dudx(0:nx+1,0:nz+1),dvdz(0:nx+1,0:nz+1),dudz(0:nx+1,0:nz+1),dvdx(0:nx+1,0:nz+1), &
            Phi(0:nx+1,0:nz+1),uold(0:nx+1,0:nz+1),vold(0:nx+1,0:nz+1),                      &
            u_previt(0:nx+1,0:nz+1),v_previt(0:nx+1,0:nz+1),                                 &
            uextr(0:nx+1,0:nz+1),vextr(0:nx+1,0:nz+1),uT(nx),uB(nx),vL(nz),vR(nz),           &
            bcT(nx),bcB(nx),bcL(nz),bcR(nz))
   
!-----------------------
! Solver allocations   
!----------------------------------------------------------------------------------   
! NB: nnz is the number of non-zero elements of the whole matrix.
! nnz => integer. Number of expected non-zero elements in the matrix
! nnzval => integer, array(nnz). Non-zero elements stored row after row
! columns => integer, array(nnz). Column index of nonzero elements
! rowind => integer, array(ndim+1). Array containing the incremental numbers 
!                                   of nonzero elements that appear in a new row
!----------------------------------------------------------------------------------
   nnz=13*(nx+1)*nz+13*(nz+1)*nx+4*(nz+1)*(nx+1)
   allocate(nnzval(nnz),columns(nnz),rowind(ndim+1), &
            irn_m(nnz),jcn_m(nnz))
   
 end subroutine stokes_alloc

 
!===========================================================================
 subroutine stokes
! Solution to the incompressible Stokes equation in 2-D Cartesian geometry
! Finite volume formulation of 'Gerya & Yuen, PEPI 140(2003)'
!
! NB Continuity equation for pressure is NOT solved at the upper left corner
!    to ensure uniqueness of the pressure field
!===========================================================================
   use share; use domain; use linsolver; use markeraux
   implicit none
   integer :: i,j,l,jmin,jmax,iu,iv,ip,innz
   real *8 :: start_time, end_time
   real (kind=8) :: diff
   real (kind=8), dimension(:,:), allocatable :: mm_matrix
   character (len=32) :: varfmt
!===========================================================================
    
   loadx=0d0 
   loadz=0d0

   nli=0
!*************************************
non_linear_iteration: do 
!*************************************

!--------------------------------------------------------
! Set strain-rate and velocity for non-linear iterations
!--------------------------------------------------------
   if(nli.gt.0) then
     epsII_previt=epsII
     u_previt=u
     v_previt=v
   endif

!------------------------
! Setup the viscosity
!------------------------
   call varparm

!-----------------------------------------------
! Main loop to setup matrix and right hand side
!-----------------------------------------------
   rhs=0d0 
   nnzval=0d0
   columns=0
   rowind=0
   innz=0

   l=0
   do i=0,nx+1

     if(i.eq.0) then
       jmin=1
       jmax=nz
     else
       jmin=0
       jmax=nz
     endif
     
     do j=jmin,jmax

!----------------------------------------------------
! Average viscosities needed for the discretization
!----------------------------------------------------
       if(i.ne.nx+1)            etaN=(eta(i,j)+eta(i,j+1)+eta(i+1,j)+eta(i+1,j+1))/4
       if(i.ne.nx+1.and.j.ne.0) etaS=(eta(i,j-1)+eta(i,j)+eta(i+1,j-1)+eta(i+1,j))/4
       if(i.ne.nx+1)            etaE=(eta(i,j)+eta(i,j+1)+eta(i+1,j)+eta(i+1,j+1))/4
       if(i.ne.0)               etaW=(eta(i-1,j)+eta(i-1,j+1)+eta(i,j)+eta(i,j+1))/4 
       
!-------------------------------------------------------------
! Load (Ra*(T-Tref)-RaC*Gamma) for non-dimensional equations
!-------------------------------------------------------------
       if(i.ne.nx+1) loadz=drho(i,j)

!------------------------------------------------------------------ 
! Momentum equation along z for v-nodes next to the left boundary
!------------------------------------------------------------------ 
       if(i.eq.0) then
         call mom_v_bcl(i,j,l,innz) 

!---------------------------------------------------------------       
! Continuity equation for p-nodes next to the right boundary
!---------------------------------------------------------------       
       elseif(i.eq.nx+1) then
         call cont_p_bcr(i,j,l,innz)

!----------------------------------------------------------------       
! Momentum and continuity equations for all other internal nodes       
!----------------------------------------------------------------       
       else
       
         if(j.eq.0) then
           call mom_u(i,j,l,innz)
           call cont_p(i,j,l,innz)
         else
           call mom_u(i,j,l,innz)
           call mom_v(i,j,l,innz)
           if(i.ne.1.or.j.ne.nz) then
           call cont_p(i,j,l,innz)
           endif
         endif

       endif
      
     end do
   end do

   rowind(ndim+1)=innz+1

!-----------------------------------------------------------
! Solution of the system
! sol: solution vector
! b: auxiliary vector to solve the system with multiple rhs
!-----------------------------------------------------------
   sol=0d0
   b(:,1)=rhs

   if((debug.eq.1)) then 
     call cpu_time(start_time)
   endif

   call mtrx_idxs(ndim,columns,rowind,irn_m,jcn_m)
   
! Retrieve and write on file the whole momentum matrix (only used for debug/development) 
   if(write_mm_matrix.eq.1) then
     print *, 'Store momentum matrix in mm_matrix.dat'      

     allocate(mm_matrix(ndim,ndim))      
     call mtrx_restore(nnzval,columns,rowind,mm_matrix,ndim)

     open(1,file='mm_matrix.dat')
     write(varfmt,99) ndim
     do i=1,ndim
        write(1,varfmt) (mm_matrix(i,j), j=1,ndim) 
     end do
     deallocate(mm_matrix)
99   format ('(',I0,'E12.3)')
   endif  

!----------------------------------------
!     Solution with MUMPS
!----------------------------------------
   mumps_par%JOB = -1
   mumps_par%SYM = 0
   mumps_par%PAR = 1
   call dmumps(mumps_par)
   mumps_par%ICNTL(1) = 0
   mumps_par%ICNTL(2) = 0
   mumps_par%ICNTL(3) = minfo
   mumps_par%ICNTL(4) = 2
   mumps_par%ICNTL(7) = mord
   mumps_par%ICNTL(23) = 4000
   mumps_par%N = ndim
   mumps_par%NZ = nnz      
   allocate(mumps_par%IRN ( mumps_par%NZ ) )
   allocate(mumps_par%JCN ( mumps_par%NZ ) )
   allocate(mumps_par%A( mumps_par%NZ ) )
   allocate(mumps_par%RHS ( mumps_par%N  ) )

   do i = 1, mumps_par%NZ
     mumps_par%IRN(i) = irn_m(i)
     mumps_par%JCN(i) = jcn_m(i)
     mumps_par%A(i) = nnzval(i)
   end do

   do i = 1, mumps_par%N
     mumps_par%RHS(i) = rhs(i)
   end do

!  Call package for solution
   mumps_par%JOB = 6
   CALL DMUMPS(mumps_par)

!  Store the mumps solution 
   sol = mumps_par%RHS

!  Deallocate user data
   deallocate( mumps_par%IRN )
   deallocate( mumps_par%JCN )
   deallocate( mumps_par%A   )
   deallocate( mumps_par%RHS )

!  Destroy the instance (deallocate internal data structures)
   mumps_par%JOB = -2
   CALL DMUMPS(mumps_par)

   if((debug.gt.1)) then 
     call cpu_time(end_time)
     write(*,'(A26,2E15.5)') '   CPU time momentum solver:', end_time-start_time
   endif

!----------------------------------------------------------
! Provisional storage of the solution on the original grid
! upom, iu: horizontal flow and index
! vpom, iv: vertical flow and index
! ppom, ip: pressure and index
!----------------------------------------------------------
!   open(1,file='solvec')
   l=0      
   upom=0d0; iu=0
   vpom=0d0; iv=0
   ppom=0d0; ip=0
   do i=0,nx+1
     if(i.eq.0) then
       jmin=1
       jmax=nz
     else
       jmin=0
       jmax=nz
     endif
     
     do j=jmin,jmax
     
       if(i.eq.0) then
         l=l+1; iv=iv+1; vpom(iv)=sol(l)         
       elseif(i.eq.nx+1) then
         l=l+1; ip=ip+1; ppom(ip)=sol(l)
       else
         if(j.eq.0) then
           l=l+1; iu=iu+1; upom(iu)=sol(l)
           l=l+1; ip=ip+1; ppom(ip)=sol(l)
         else
           l=l+1; iu=iu+1; upom(iu)=sol(l)
           l=l+1; iv=iv+1; vpom(iv)=sol(l)
           if(i.ne.1.or.j.ne.nz) then
           l=l+1; ip=ip+1; ppom(ip)=sol(l)
           endif
         endif
       endif
       
     end do
   end do
!  close(1)

!--------------------------------------------------------------------
! From the solution vector 'sol' write both components of the flow 
! on the original grid and pressure on its own grid along with the
! vertical stress
!--------------------------------------------------------------------
   call writesol

!--------------------------------------------------------
! For linear viscosity, exit the iteration loop
!--------------------------------------------------------
   if(nlv.eq.0) exit

!-------------------------------------------------------------------------------------
! For the first non-linear iteration, increment the iteration index without computing 
! the difference with respect to the previous solution and cycle
!-------------------------------------------------------------------------------------
   if(nli.eq.0) then
     nli=nli+1
     cycle
   endif

!-----------------------------------------------------------------------------------
! For the subsequent non-linear iterations, compute the difference and check
! convergence criterion
!-----------------------------------------------------------------------------------
! based on velocity
   if(conv_test.eq.1) then
     diff=max(maxval(dabs(u-u_previt)/maxval(u_previt)),maxval(dabs(v-v_previt)/maxval(v_previt)))
! or based on strain-rate
   elseif(conv_test.eq.2) then
     diff=maxval(dabs(epsII-epsII_previt)/maxval(epsII_previt))
   endif

   if((debug.eq.1).and.(nlv.ne.0)) write(*,'(A18,I3,A6,E14.5)') 'Stokes iteration:',nli,'diff=',diff
   if(((nlv.ne.0).and.(nli.ge.1).and.(diff.le.nlitol)).or.(nli.gt.maxit)) exit

   nli=nli+1

!******************************* 
   end do non_linear_iteration
!******************************* 

   eta_old=eta

 end subroutine stokes  


!=====================================================================
 subroutine varparm
! Setup variable parameters on the grid: 
! viscosity, conductivity, expansivity, RHS of Stokes equation
!=====================================================================
   use share; use domain; use markeraux
   implicit none
   integer :: i,j,k
   real (kind=8) :: x,z,T,f,e,ediff,eyield,edisl,z0,rho_r,drho_r,T_r,zr, &
                    fph,etac,Tph,Cl,arr,phase_bu,n,ae_ref,av_ref,eps
   real (kind=8), dimension(3) :: vm    
   real (kind=8), dimension(:), allocatable :: ff    
   fph(z)=0.5d0*(1d0+dtanh((z-z0-Cl*(T-Tph))/dph))
!=====================================================================

    if(nphases.ne.0) allocate(ff(nphases))

!---------------------------------------------
! Setup the viscosity, thermal conductivity
! and thermal expansivity
!---------------------------------------------
   do i=0,nx+1 
     x=xgrid(i)

     do j=0,nz+1
       z=dimz-zgrid(j)
       T=Tnew(i,j)

!----------------------------------------------
!                Viscosity
!----------------------------------------------

! Initialize array of viscosity mechanism to a very large value that will be
! overwritten from the viscosity associated with the actual deformation mechanisms at work
       vm=1d100

!********************************************************************************
! Compute phase functions and, if required, set different activation parameters
! for each phase
!********************************************************************************
       if(nphases.ne.0) then
         do k=1,nphases
           z0=phase_z(k)
           Cl=phase_Cl(k)
           Tph=phase_T(k)
           ff(k)=fph(z)
         end do

         if(phase_dep_ap.eq.1) then

           if(visc_law.eq.1) then
             print *, "Phase-dependent activation parameters cannot be used yet with FK viscosity law. Use Arrhenius instead."
             print *, "Cheers"
             stop
           endif

           do k=1,nphases

             if (k.eq.1.and.ff(k).le.0.5d0) then 
               ae = phase_ae(k)
               av = phase_av(k)
               ae_ds = phase_ae_ds(k)
               av_ds = phase_av_ds(k)
               prefac = phase_prefac(k)
               prefac_ds = phase_prefac_ds(k)
               exit
             endif  

             if (k.eq.nphases.and.ff(k).gt.0.5d0) then 
               ae = phase_ae(k+1)
               av = phase_av(k+1)
               ae_ds = phase_ae_ds(k+1)
               av_ds = phase_av_ds(k+1)
               prefac = phase_prefac(k+1)
               prefac_ds = phase_prefac_ds(k+1)
               exit
             endif  

             if (ff(k).gt.0.5d0.and.ff(k+1).le.0.5d0) then 
               ae = phase_ae(k+1)
               av = phase_av(k+1)
               ae_ds = phase_ae_ds(k+1)
               av_ds = phase_av_ds(k+1)
               prefac = phase_prefac(k+1)
               prefac_ds = phase_prefac_ds(k+1)
               exit
             endif

           end do  


         endif

       endif

       if(nphases.gt.0) then
         ae_ref = phase_ae(1)
         av_ref = phase_av(1)
       else
         ae_ref = ae
         av_ref = av
       endif

!***********************************
! T- and p- dependence
!***********************************
! Frank-Kamentzkii law
       if(visc_law.eq.1) then
!       Temperature dependent part
         f=dexp(dlog(Ft)*(Tref-Tnew(i,j)))
       
!       Pressure dependent part
         ediff=dexp(dlog(Fz)*(z-zref))*f
         vm(1)=ediff 
         
! Arrhenius law for diffusion creep
       elseif(visc_law.eq.2) then
         arr=(ae_ref + av_ref*zref) / (Tref + Tsurf)
         ediff=prefac * exp((ae + av*z) / (Tnew(i,j) + Tsurf) - arr)
         vm(1)=ediff 

! Viscosity calculated from particles in mkadvect
       elseif(visc_law.eq.4) then
         ediff=eta(i,j)
         vm(1)=ediff 

       endif

!*****************************************************
! Non-linear viscosity (mixed or purely non-linear)
!*****************************************************
       if(nlv.ne.0) then

!      To compute the viscosity for first non linear iteration and first time
!      step, a reference value of the strain rate is used
         if(nli.eq.0.and.ti.eq.0) then
           eps=epsref
         else
           eps=epsII(i,j)
         endif

!      Plastic viscosity               
         if(((sigmay.ne.0d0).or.(sigmayz.ne.0d0)).and.z.le.yield_depth) then
           eyield=(sigmay + sigmayz*z)/(2d0*eps)
!           eyield=1d-3+1d0/eps
         else
           eyield=1d100
         endif
         vm(2)=eyield 

!      Dislocation creep viscosity
         if(stress_exp.gt.0d0) then
           n=stress_exp

           if(visc_law.eq.2) then
!            non-dimensional dislocation-creep viscosity for mixed rheology (normalized by diffusion creep viscosity)
             arr=(ae_ref + av_ref*zref) / (Tref + Tsurf)
             edisl=prefac_ds * exp((ae_ds + av_ds*z) / (n*(Tnew(i,j) + Tsurf)) - arr) * eps**((1d0-n)/n)     

           elseif(visc_law.eq.3) then  
!            non-dimensional dislocation-creep viscosity for purely dislocation models
             arr=(ae_ds + av_ds*zref) / (n*(Tref + Tsurf))
             edisl=epsref**((n-1d0)/n) * exp((ae_ds + av_ds*z) / (n*(Tnew(i,j) + Tsurf)) - arr) * eps**((1d0-n)/n)
           endif  
         else
           edisl=1d100
         endif
         vm(3)=edisl

!       The first non-linear iteration uses the viscosity from the previous time-step
         if(nli.eq.0.and.ti.gt.0) then
           e=eta_old(i,j)

         else

!         Mixed rheology (FK diffusion + plasticity)
           if(visc_law.eq.1) then
             e=1d0/(1d0/eyield + 1d0/ediff)

!         Mixed rheology (Arrhenius diffusion + dislocation + plasticity)
           elseif(visc_law.eq.2) then
             e=1d0/(1d0/eyield + 1d0/edisl + 1d0/ediff)

!         Purely non-linear rheology (dislocation + plasticity)
           elseif(visc_law.eq.3) then  
             e=1d0/(1d0/eyield + 1d0/edisl)

           endif

         endif

!********************
! Linear viscosity
!********************
       else
         e=ediff

       endif

!***********************************
! Phase dependent jumps
!***********************************
       if(nphases.ne.0) then

         do k=1,nphases
           etac=phase_eta(k)
           e=e+ff(k)*(e*etac-e)
         end do

       endif

! Set the viscosity on the grid      
       eta(i,j)=e

! Determine the dominant viscosity mechanism at grid point i,j and store it as integer into visc_mech
       visc_mech(i,j)=minval(minloc(vm))

!----------------------------------------------
!           Thermal expansivity       
!----------------------------------------------
! Temperature-dependent alpha
       if(alpha_temp.eq.1) then
         T=Tnew(i,j)
! No temperature dependence: the formula alpha(p,T) is evaluated
! along a linear temperature profile
       elseif(alpha_temp.eq.0) then
         T=z
       endif

! alpha_law=0 -> constant 
       if(alpha_law.eq.0) then 
         alpha(i,j)=1d0

! alpha_law=1 -> pT-dependent alpha (same parameters for all phases)
       elseif(alpha_law.eq.1) then 
         call refstate(zgrid(j),rho_r,drho_r,T_r)

         alpha(i,j) = (alpha0 + alpha1*(Tsurf+T) + alpha2/(Tsurf+T)**2d0) * exp(-alpha3*z)
!         alpha(i,j) = alpha(i,j)/(rho_r)**alpha3

! alpha_law=2 -> pT-dependent alpha (with variable coefficients for each phase)
       elseif(nphases.ne.0.and.alpha_law.eq.2) then 

         do k=1,nphases

           if (k.eq.1.and.ff(k).le.0.5d0) then 
             alpha(i,j) = (phase_alpha0(k) + phase_alpha1(k)*(Tsurf+T) + phase_alpha2(k)/(Tsurf+T)**2d0) &
                          *exp(-phase_alpha3(k)*z)
             exit
           endif  

           if (k.eq.nphases.and.ff(k).gt.0.5d0) then 
             alpha(i,j) = (phase_alpha0(k+1) + phase_alpha1(k+1)*(Tsurf+T) + phase_alpha2(k+1)/(Tsurf+T)**2d0) &
                          *exp(-phase_alpha3(k+1)*z)
             exit
           endif  

           if (ff(k).gt.0.5d0.and.ff(k+1).le.0.5d0) then 
             alpha(i,j) = (phase_alpha0(k+1) + phase_alpha1(k+1)*(Tsurf+T) + phase_alpha2(k+1)/(Tsurf+T)**2d0) &
                          *exp(-phase_alpha3(k+1)*z)
             exit
           endif

         end do  

       endif

!----------------------------------------------
!          Thermal conductivity
!----------------------------------------------
! Temperature-dependent k
       if(kappa_temp.eq.1) then
         T=Tnew(i,j)
! No temperature dependence: the formula for k(p,T) is evaluated
! along a linear temperature profile
       elseif(kappa_temp.eq.0) then
         T=z
       endif

! kappa_law=0 -> Constant thermal conductivity
       if(kappa_law.eq.0) then 
         kap(i,j)=1d0

! kappa_law=1 -> pT-dependent conductivity (same parameters for all phases) 
       elseif(kappa_law.eq.1) then 
         kap(i,j) = (kap0 + kap1*z)*(1d0+T/Tsurf)**kap2

! kappa_law=2 -> pT-dependent alpha (with variable coefficients for each phase)
       elseif(nphases.ne.0.and.kappa_law.eq.2) then 

         do k=1,nphases

           if (k.eq.1.and.ff(k).le.0.5d0) then 
             kap(i,j) = (phase_kap0(k) + phase_kap1(k)*z)*(1d0+T/Tsurf)**phase_kap2(k)
             exit
           endif

           if (k.eq.nphases.and.ff(k).gt.0.5d0) then 
             kap(i,j) = (phase_kap0(k+1) + phase_kap1(k+1)*z)*(1d0+T/Tsurf)**phase_kap2(k+1)
             exit
           endif

           if (ff(k).gt.0.5d0.and.ff(k+1).le.0.5d0) then 
             kap(i,j) = (phase_kap0(k+1) + phase_kap1(k+1)*z)*(1d0+T/Tsurf)**phase_kap2(k+1)
             exit
           endif


         end do
       
       endif 

     end do  
   end do

! Overwrite part of the viscosity field if required
   if(overwr_visc.eq.1) then
      call overwrite_viscosity
   endif

! Limit viscosity range
   do i=0,nx+1 
     do j=0,nz+1
       if(eta(i,j).lt.etamin) then
         eta(i,j)=etamin
         visc_mech(i,j)=4
       endif

       if(eta(i,j).gt.etamax) then
         eta(i,j)=etamax
         visc_mech(i,j)=5
       endif
     end do
   end do

!------------------------------------------
! Setup the buoyancy term
!------------------------------------------
   do i=0,nx+1 
     do j=0,nz+1
       z=dimz-zgrid(j)

       if(cpr.eq.1) then
         zr=zgrid(j)
         call refstate(zr,rho_r,drho_r,T_r)
       else
         rho_r=1d0
         drho_r=1d0
         T_r=0d0
       endif

!*********************************
! Compute phase functions 
!*********************************
       phase_bu=0d0
       if(nphases.ne.0) then
         T=Tnew(i,j)
         
         do k=1,nphases
           z0=phase_z(k)
           Cl=phase_Cl(k)
           Tph=phase_T(k)
           ff(k)=fph(z)
           phase_bu=phase_bu+ff(k)*phase_Rb(k)
           phase_f(k,i,j)=ff(k)
         end do

       endif

       if(i.ne.nx+1) drho(i,j)=-rho_r*alpha(i,j)*((Tnew(i,j)+Tnew(i+1,j))/2d0-T_r)*Ra  &
                               +(Gam(i,j)+Gam(i+1,j))/2d0*RaC                          &
                               +phase_bu

     end do
   end do

   if(nphases.ne.0) deallocate(ff)

       
 end subroutine varparm

!============================================================================
 subroutine overwrite_viscosity
! Overwrite (part of) the viscosity field obtained from subroutine varparm  
! with a predefined viscosity distribution
!============================================================================
   use share; use domain; use markeraux
   implicit none
   integer :: i,j
   real (kind=8) :: x,z,zslow,zsup,x1,z1,x2,z2,x3,z3,x4,z4,x_trench
   real (kind=8), parameter :: eta_seed=1d0, seed_width=0.027d0, seed_depth=0.027d0, dip_angle=pi/6d0
   character (len=100) :: readval
!=======================================================================
 
   call readinput("x_trench",readval)
   read(readval,*) x_trench

   if(overwr_visc_mod.eq.1) then
     x1=x_trench; z1=0d0
     x2=x_trench+seed_width; z2=0d0
     x3=x2+seed_depth*dtan(pi/2d0-dip_angle); z3=seed_depth
     x4=x1+seed_depth*dtan(pi/2d0-dip_angle); z4=seed_depth
     
     do i=0,nx+1
       x=xgrid(i)

       zslow=z1+(x1-x)/(x1-x4)*(z4-z1)
       zsup=z2+(x2-x)/(x2-x3)*(z3-z2)
       
       do j=0,nz+1
         z=zgrid(j)
 
         if(((dimz-z).le.seed_depth).and.((dimz-z).le.zslow).and.((dimz-z).ge.zsup)) then
           eta(i,j)=eta_seed 
         endif

       end do
     end do

   else
     print *, "overwrite_viscosity_model=",overwr_visc_mod,"does not exist"
     print *, "Cheers"
     stop

   endif  


 end subroutine overwrite_viscosity


!==============================================================
 subroutine mom_v_bcl(i,j,l,innz)
! Momentum equation along z for nodes v(i+1/2,j) at the right 
! of the left boundary    4
! NB the routine is always accessed when i=0 (and j=1,...,nz)
!==============================================================
   use share; use domain
   implicit none
   integer, intent(in) :: i,j
   integer, intent(inout) :: l,innz
   real (kind=8) :: z,rho_r,drho_r,T_r
!===========================================================
   
   if(cpr.eq.1) then
     z=zgrid(j)+delz(j+1)/2
     call refstate(z,rho_r,drho_r,T_r)
   else
     rho_r=1d0
     drho_r=1d0
     T_r=0d0
   endif
  
   l=l+1
   rowind(l)=innz+1
   
   if(j.eq.1) then

     innz=innz+1
     nnzval(innz)=-2/(delz(j)+delz(j+1))*(2*etaN/delz(j+1)+2*etaS/delz(j))         &  !v(i+1/2,j)
                  -2*eta(i+1,j)/delx(i+1)/(delx(i+1)+delx(i+2))                    &
                  -bcL(j)*2*eta(i,j)*2/delx(i+1)/(delx(i+1)+delx(i+1))             &
                  +cpr*2/(delz(j)+delz(j+1))*2d0/3d0*(etaN/delz(j+1)+etaS/delz(j))
     columns(innz)=l
     
     innz=innz+1
     nnzval(innz)=2/(delz(j)+delz(j+1))*2*etaN/delz(j+1)            &                 !v(i+1/2,j+1)
                  -cpr*2/(delz(j)+delz(j+1))*2d0/3d0*etaN/delz(j+1)
     columns(innz)=l+1
     
     innz=innz+1
     nnzval(innz)=-1/delx(i+1)*2*eta(i+1,j)/(delz(j)+delz(j+1))     &                 !u(i+1,j-1/2)
                  +cpr*2/(delz(j)+delz(j+1))*2d0/3d0*etaS/delx(i+1)
     columns(innz)=l+nz
     
     innz=innz+1
     nnzval(innz)=2/(delz(j)+delz(j+1))-ala*Di*rho_r/2d0                              !p(i+1/2,j-1/2)
     columns(innz)=l+(nz+1)
     
     innz=innz+1
     nnzval(innz)=1/delx(i+1)*2*eta(i+1,j)/(delz(j)+delz(j+1))      &                 !u(i+1,j+1/2)
                  -cpr*2/(delz(j)+delz(j+1))*2d0/3d0*etaN/delx(i+1)
     columns(innz)=l+(nz+2)
     
     innz=innz+1
     nnzval(innz)=1/delx(i+1)*2*eta(i+1,j)/(delx(i+1)+delx(i+2))                      !v(i+3/2,j)
     columns(innz)=l+(nz+3)
     
     innz=innz+1
     nnzval(innz)=-2/(delz(j)+delz(j+1))-ala*Di*rho_r/2d0                             !p(i+1/2,j+1/2)
     columns(innz)=l+(nz+4)
     
   else  
     
     innz=innz+1
     nnzval(innz)=2/(delz(j)+delz(j+1))*2*etaS/delz(j)            &                   !v(i+1/2,j-1) 
                  -cpr*2/(delz(j)+delz(j+1))*2d0/3d0*etaS/delz(j)
     columns(innz)=l-1
   
     innz=innz+1
     nnzval(innz)=-2/(delz(j)+delz(j+1))*(2*etaN/delz(j+1)+2*etaS/delz(j))          & !v(i+1/2,j)
                  -1/delx(i+1)*2*eta(i+1,j)/(delx(i+1)+delx(i+2))                   &
                  -bcL(j)*2*eta(i,j)*2/delx(i+1)/(delx(i+1)+delx(i+1))              &
                  +cpr*2/(delz(j)+delz(j+1))*2d0/3d0*(etaN/delz(j+1)+etaS/delz(j))
     columns(innz)=l
     
     if(j.lt.nz) then
       innz=innz+1
       nnzval(innz)=2/(delz(j)+delz(j+1))*2*etaN/delz(j+1)            &               !v(i+1/2,j+1)
                    -cpr*2/(delz(j)+delz(j+1))*2d0/3d0*etaN/delz(j+1)
       columns(innz)=l+1
     endif
     
     innz=innz+1
     nnzval(innz)=-1/delx(i+1)*2*eta(i+1,j)/(delz(j)+delz(j+1))      &                !u(i+1,j-1/2) 
                  +cpr*2/(delz(j)+delz(j+1))*2d0/3d0*etaS/delx(i+1)
     columns(innz)=l+(nz+2*j-3)
     
     innz=innz+1
     nnzval(innz)=2/(delz(j)+delz(j+1))-ala*Di*rho_r/2d0                              !p(i+1/2,j-1/2)
     columns(innz)=l+(nz+2*j-1)
     
     innz=innz+1
     nnzval(innz)=1/delx(i+1)*2*eta(i+1,j)/(delz(j)+delz(j+1))      &                 !u(i+1,j+1/2)
                  -cpr*2/(delz(j)+delz(j+1))*2d0/3d0*etaN/delx(i+1)
     columns(innz)=l+(nz+2*j)
     
     innz=innz+1
     nnzval(innz)=1/delx(i+1)*2*eta(i+1,j)/(delx(i+1)+delx(i+2))                      !v(i+3/2,j)
     columns(innz)=l+(nz+2*j+1)
     
     if(j.ne.nz) then
       innz=innz+1
       nnzval(innz)=-2/(delz(j)+delz(j+1))-ala*Di*rho_r/2d0                           !p(i+1/2,j+1/2)
       columns(innz)=l+(nz+2*j+2)
     endif

   endif      

   rhs(l)=loadz+bcL(j)*(-1/delx(i+1)*2*eta(i,j)*2/(delx(i+1)+delx(i+1))*vL(j))


 end subroutine mom_v_bcl
  
!===========================================================
 subroutine cont_p_bcr(i,j,l,innz)
! Continuity equation for nodes p(i-1/2,j-1/2) at the left 
! of the right boundary
! NB The subroutine is always accessed with i=nx+1
!===========================================================
   use share; use domain
   implicit none
   integer, intent(in) :: i,j
   integer, intent(inout) :: l,innz
   real (kind=8) :: z,rho_r,drho_r,T_r
!===========================================================
   
   if(cpr.eq.1) then
     z=zgrid(j)+delz(j+1)/2
     call refstate(z,rho_r,drho_r,T_r)
   else
     rho_r=1d0
     drho_r=1d0
     T_r=0d0
   endif
  
   l=l+1
   rowind(l)=innz+1

!   if(rgx.eq.0.or.rgz.eq.0) then

   if(j.eq.0) then

     innz=innz+1
     nnzval(innz)=-rho_r*1/delx(i)         !u(i-1,j-1/2) 
     columns(innz)=l-(3*nz+2)

     innz=innz+1
     nnzval(innz)=rho_r*1/delz(j+1)+   &   !v(i-1/2,j)
                  cpr*drho_r/2  
     columns(innz)=l-(3*nz-1)

   else

     innz=innz+1
     nnzval(innz)=-rho_r*1/delx(i)         !u(i-1,j-1/2) 
     columns(innz)=l-(3*nz+3-2*j)

     innz=innz+1
     nnzval(innz)=-rho_r*1/delz(j+1)+  &   !v(i-1/2,j-1) 
                  cpr*drho_r/2  
     columns(innz)=l-(3*nz+2-2*j)

     if(j.ne.nz) then
     innz=innz+1
     nnzval(innz)=rho_r*1/delz(j+1)+   &   !v(i-1/2,j)     
                  cpr*drho_r/2  
     columns(innz)=l-(3*nz-1-2*j)
     endif

   endif

!   endif

!-------------------------------------------------------------
!  NB: The row built in this subroutine has zero element 
!      on the diagonal that must be stored anyway 
!-------------------------------------------------------------
   innz=innz+1
   nnzval(innz)=0d0
   columns(innz)=l

 end subroutine cont_p_bcr


!==================================================
 subroutine mom_u(i,j,l,innz)
! Momentum equation along x for nodes u(i,j+1/2)
!==================================================
   use share; use domain
   implicit none
   integer :: k,kS
   integer, intent(in) :: i,j
   integer, intent(inout) :: l,innz
!==================================================

   k=1
   if(j.eq.nz) k=0
   kS=0
   if(j.eq.nz) kS=1


   l=l+1
   rowind(l)=innz+1

   if(j.eq.0.and.i.eq.1) then  
 
     innz=innz+1
     nnzval(innz)=-1/delz(j+1)*2*eta(i,j+1)/(delx(i)+delx(i+1))     &              !v(i-1/2,j+1)
                  +cpr*2/(delx(i)+delx(i+1))*2d0/3d0*etaW/delz(j+1) 
     columns(innz)=l-nz
!     endif

     innz=innz+1
     nnzval(innz)=-2/(delx(i)+delx(i+1))*(2*etaE/delx(i+1)+2*etaW/delx(i)) &       !u(i,j+1/2)
                  -1/delz(j+1)*2*eta(i,j+1)/(delz(j+1)+delz(j+2))          &
                  -1/delz(j+1)*bcB(i)*2*eta(i,j)*2/(delz(j+1)+delz(j+1))      &
                  +cpr*2/(delx(i)+delx(i+1))*2d0/3d0*(etaE/delx(i+1)+etaW/delx(i))
     columns(innz)=l
     
     innz=innz+1
     nnzval(innz)=2/(delx(i)+delx(i+1))                                   !p(i-1/2,j+1/2) 
     columns(innz)=l+1
     
     innz=innz+1
     nnzval(innz)=1/delz(j+1)*2*eta(i,j+1)/(delz(j+1)+delz(j+2))          !u(i,j+3/2) 
     columns(innz)=l+2
     
     innz=innz+1
     nnzval(innz)=1/delz(j+1)*2*eta(i,j+1)/(delx(i)+delx(i+1))      &     !v(i+1/2,j+1) 
                  -cpr*2/(delx(i)+delx(i+1))*2d0/3d0*etaE/delz(j+1) 
     columns(innz)=l+3
     
     innz=innz+1
     nnzval(innz)=2/(delx(i)+delx(i+1))*2*etaE/delx(i+1)           &       !u(i+1,j+1/2)
                 -cpr*2/(delx(i)+delx(i+1))*2d0/3d0*etaE/delx(i+1)
     columns(innz)=l+(3*nz+1)
     
     innz=innz+1
     nnzval(innz)=-2/(delx(i)+delx(i+1))                                   !p(i+1/2,j+1/2) 
     columns(innz)=l+(3*nz+2)

     rhs(l)=loadx+bcB(i)*(-1/delz(j+1)*2*eta(i,j)*2/(delz(j+1)+delz(j+1))*uB(i))

   elseif(j.eq.0.and.i.ne.1) then  
 
     if(i.eq.2) then
       innz=innz+1
       nnzval(innz)=2/(delx(i)+delx(i+1))*2*etaW/delx(i)            &      !u(i-1,j+1/2) *
                    -cpr*2/(delx(i)+delx(i+1))*2d0/3d0*etaW/delx(i)
       columns(innz)=l-(3*nz+1)
     else
       innz=innz+1
       nnzval(innz)=2/(delx(i)+delx(i+1))*2*etaW/delx(i)            &      !u(i-1,j+1/2) *
                    -cpr*2/(delx(i)+delx(i+1))*2d0/3d0*etaW/delx(i)
       columns(innz)=l-(3*nz+2)
     endif

     if(i.eq.2) then
       innz=innz+1
       nnzval(innz)=-1/delz(j+1)*2*eta(i,j+1)/(delx(i)+delx(i+1))     &    !v(i-1/2,j+1) !
                    +cpr*2/(delx(i)+delx(i+1))*2d0/3d0*etaW/delz(j+1) 
       columns(innz)=l-(3*nz-2)
     else
       innz=innz+1
       nnzval(innz)=-1/delz(j+1)*2*eta(i,j+1)/(delx(i)+delx(i+1))     &    !v(i-1/2,j+1) !
                    +cpr*2/(delx(i)+delx(i+1))*2d0/3d0*etaW/delz(j+1) 
       columns(innz)=l-(3*nz-1)
     endif

     innz=innz+1
     nnzval(innz)=-2/(delx(i)+delx(i+1))*(2*etaE/delx(i+1)+2*etaW/delx(i)) & !u(i,j+1/2)
                  -1/delz(j+1)*2*eta(i,j+1)/(delz(j+1)+delz(j+2))          &
                  -1/delz(j+1)*bcB(i)*2*eta(i,j)*2/(delz(j+1)+delz(j+1))      &
                  +cpr*2/(delx(i)+delx(i+1))*2d0/3d0*(etaE/delx(i+1)+etaW/delx(i))
     columns(innz)=l
     
     innz=innz+1
     nnzval(innz)=2/(delx(i)+delx(i+1))                                !p(i-1/2,j+1/2) 
     columns(innz)=l+1
     
     innz=innz+1
     nnzval(innz)=1/delz(j+1)*2*eta(i,j+1)/(delz(j+1)+delz(j+2))       !u(i,j+3/2) 
     columns(innz)=l+2
     
     innz=innz+1
     nnzval(innz)=1/delz(j+1)*2*eta(i,j+1)/(delx(i)+delx(i+1))       & !v(i+1/2,j+1) 
                  -cpr*2/(delx(i)+delx(i+1))*2d0/3d0*etaE/delz(j+1) 
     columns(innz)=l+3
     
     if(i.ne.nx) then
       innz=innz+1
       nnzval(innz)=2/(delx(i)+delx(i+1))*2*etaE/delx(i+1)            & !u(i+1,j+1/2)
                    -cpr*2/(delx(i)+delx(i+1))*2d0/3d0*etaE/delx(i+1)
       columns(innz)=l+(3*nz+2)
     endif
     
     if(i.ne.nx) then
       innz=innz+1
       nnzval(innz)=-2/(delx(i)+delx(i+1))                         !p(i+1/2,j+1/2) *
       columns(innz)=l+(3*nz+3)
     endif
     
     if(i.eq.nx) then
       innz=innz+1 
       nnzval(innz)=-2/(delx(i)+delx(i+1))                         !p(i+1/2,j+1/2) *
       columns(innz)=l+(3*nz+2)
     endif

     rhs(l)=loadx+bcB(i)*(-1/delz(j+1)*2*eta(i,j)*2/(delz(j+1)+delz(j+1))*uB(i))
    
   else

     if(i.eq.1) then

       innz=innz+1
       nnzval(innz)=1/delz(j+1)*2*eta(i,j)/(delx(i)+delx(i+1))       &  !v(i-1/2,j)   *
                    -cpr*2/(delx(i)+delx(i+1))*2d0/3d0*etaW/delz(j+1) 
       columns(innz)=l-(nz+2*j)

       if(j.ne.nz) then
       innz=innz+1
       nnzval(innz)=-1/delz(j+1)*2*eta(i,j+1)/(delx(i)+delx(i+1))     & !v(i-1/2,j+1) +
                    +cpr*2/(delx(i)+delx(i+1))*2d0/3d0*etaW/delz(j+1) 
       columns(innz)=l-(nz+2*j-1)
       endif

     elseif(i.eq.2) then

       innz=innz+1
       nnzval(innz)=2/(delx(i)+delx(i+1))*2*etaW/delx(i)            &   !u(i-1,j+1/2) -
                    -cpr*2/(delx(i)+delx(i+1))*2d0/3d0*etaW/delx(i)
       columns(innz)=l-(3*nz+1)

       innz=innz+1
       nnzval(innz)=1/delz(j+1)*2*eta(i,j)/(delx(i)+delx(i+1))        & !v(i-1/2,j)   *
                    -cpr*2/(delx(i)+delx(i+1))*2d0/3d0*etaW/delz(j+1) 
       columns(innz)=l-(3*nz)

       if(j.ne.nz) then
       innz=innz+1
       nnzval(innz)=-1/delz(j+1)*2*eta(i,j+1)/(delx(i)+delx(i+1))     & !v(i-1/2,j+1) +
                    +cpr*2/(delx(i)+delx(i+1))*2d0/3d0*etaW/delz(j+1) 
       columns(innz)=l-(3*nz-3)
       endif

     else

       innz=innz+1
       nnzval(innz)=2/(delx(i)+delx(i+1))*2*etaW/delx(i)            &   !u(i-1,j+1/2) -
                    -cpr*2/(delx(i)+delx(i+1))*2d0/3d0*etaW/delx(i)
       columns(innz)=l-(3*nz+2)

       innz=innz+1
       nnzval(innz)=1/delz(j+1)*2*eta(i,j)/(delx(i)+delx(i+1))        & !v(i-1/2,j)   *
                    -cpr*2/(delx(i)+delx(i+1))*2d0/3d0*etaW/delz(j+1) 
       columns(innz)=l-(3*nz+1)

       if(j.ne.nz) then
       innz=innz+1
       nnzval(innz)=-1/delz(j+1)*2*eta(i,j+1)/(delx(i)+delx(i+1))     &  !v(i-1/2,j+1) +
                    +cpr*2/(delx(i)+delx(i+1))*2d0/3d0*etaW/delz(j+1) 
       columns(innz)=l-(3*nz-2)
       endif

     endif

     if(j.eq.1) then 
       innz=innz+1
       nnzval(innz)=1/delz(j+1)*2*eta(i,j)/(delz(j)+delz(j+1))  !u(i,j-1/2) *
       columns(innz)=l-2
     else  
       innz=innz+1
       nnzval(innz)=1/delz(j+1)*2*eta(i,j)/(delz(j)+delz(j+1))  !u(i,j-1/2) *
       columns(innz)=l-3
     endif

     innz=innz+1
     if(j.ne.nz) then
     nnzval(innz)=-2/(delx(i)+delx(i+1))*(2*etaE/delx(i+1)+2*etaW/delx(i)) &  !u(i,j+1/2)
                  -1/delz(j+1)*2*eta(i,j+1)/(delz(j+1)+delz(j+2))          &
                  -1/delz(j+1)*2*eta(i,j)/(delz(j)+delz(j+1))              &
                  +cpr*2/(delx(i)+delx(i+1))*2d0/3d0*(etaE/delx(i+1)+etaW/delx(i))
     else
     nnzval(innz)=-2/(delx(i)+delx(i+1))*(2*etaE/delx(i+1)+2*etaW/delx(i)) &  !u(i,j+1/2)
                  -1/delz(j+1)*2*eta(i,j)/(delz(j)+delz(j+1))              &
                  -1/delz(j+1)*bcT(i)*2*eta(i,j+1)*2/(delz(j+1)+delz(j+1))    &
                  +cpr*2/(delx(i)+delx(i+1))*2d0/3d0*(etaE/delx(i+1)+etaW/delx(i))
     endif
     columns(innz)=l
     
     innz=innz+1
     nnzval(innz)=-1/delz(j+1)*2*eta(i,j)/(delx(i)+delx(i+1))       &  !v(i+1/2,j)
                  +cpr*2/(delx(i)+delx(i+1))*2d0/3d0*etaE/delz(j+1) 
     columns(innz)=l+1

     if(i.ne.1.or.j.ne.nz) then
       innz=innz+1
       nnzval(innz)=2/(delx(i)+delx(i+1))                            !p(i-1/2,j+1/2)
       columns(innz)=l+2
     endif
     
     if(j.ne.nz) then
       innz=innz+1
       nnzval(innz)=1/delz(j+1)*2*eta(i,j+1)/(delz(j+1)+delz(j+2))       !u(i,j+3/2)
       columns(innz)=l+3
       
       innz=innz+1
       nnzval(innz)=1/delz(j+1)*2*eta(i,j+1)/(delx(i)+delx(i+1))      &  !v(i+1/2,j+1)
                    -cpr*2/(delx(i)+delx(i+1))*2d0/3d0*etaE/delz(j+1) 
       columns(innz)=l+4
       
     endif

     if(i.ne.nx) then      
       if(i.eq.1) then
         innz=innz+1
         nnzval(innz)=2/(delx(i)+delx(i+1))*2*etaE/delx(i+1)            &  !u(i+1,j+1/2) +
                      -cpr*2/(delx(i)+delx(i+1))*2d0/3d0*etaE/delx(i+1)
         columns(innz)=l+(3*nz+1)
       endif
       
       if(i.ne.1) then
         innz=innz+1
         nnzval(innz)=2/(delx(i)+delx(i+1))*2*etaE/delx(i+1)            &  !u(i+1,j+1/2) +
                      -cpr*2/(delx(i)+delx(i+1))*2d0/3d0*etaE/delx(i+1)
         columns(innz)=l+(3*nz+2)
       endif
       
       if(i.eq.1) then
         innz=innz+1
         nnzval(innz)=-2/(delx(i)+delx(i+1))     !p(i+1/2,j+1/2) *
         columns(innz)=l+(3*nz+3)
       endif
       
       if(i.ne.1) then
         innz=innz+1
         nnzval(innz)=-2/(delx(i)+delx(i+1))     !p(i+1/2,j+1/2) *
         columns(innz)=l+(3*nz+4)
       endif
       
     endif
     
     if(i.eq.nx) then
       innz=innz+1
       nnzval(innz)=-2/(delx(i)+delx(i+1))       !p(i+1/2,j+1/2) *
       columns(innz)=l+(3*nz+3-2*j)
     endif
     
     rhs(l)=loadx+kS*bcT(i)*(-1/delz(j+1)*2*eta(i,j+1)*2/(delz(j+1)+delz(j+1))*uT(i))

   endif

 end subroutine mom_u

!===================================================
 subroutine mom_v(i,j,l,innz)
! Momentum equation along z for nodes v(i+1/2,j)
!===================================================
   use share; use domain
   implicit none
   integer :: k,kE
   integer, intent(in) :: i,j
   integer, intent(inout) :: l,innz
   real (kind=8) :: z,rho_r,drho_r,T_r
!===========================================================
   
   if(cpr.eq.1) then
     z=zgrid(j)+delz(j+1)/2
     call refstate(z,rho_r,drho_r,T_r)
   else
     rho_r=1d0
     drho_r=1d0
     T_r=0d0
   endif

   k=1
   if(i.eq.nx) k=0
   kE=0
   if(i.eq.nx) kE=1

   l=l+1
   rowind(l)=innz+1

   if(i.eq.1) then
     innz=innz+1
     nnzval(innz)=1/delx(i+1)*2*eta(i,j)/(delx(i)+delx(i+1)) !v(i-1/2,j) * 
     columns(innz)=l-(nz+2*j+1)
   elseif(i.eq.2) then 
     innz=innz+1
     nnzval(innz)=1/delx(i+1)*2*eta(i,j)/(delx(i)+delx(i+1)) !v(i-1/2,j) *
     columns(innz)=l-(3*nz+1)
   else 
     innz=innz+1
     nnzval(innz)=1/delx(i+1)*2*eta(i,j)/(delx(i)+delx(i+1)) !v(i-1/2,j) *
     columns(innz)=l-(3*nz+2)
   endif   

   if(j.eq.1) then 
     innz=innz+1
     nnzval(innz)=1/delx(i+1)*2*eta(i,j)/(delz(j)+delz(j+1))        & !u(i,j-1/2) = 
                  -cpr*2/(delz(j)+delz(j+1))*2d0/3d0*etaS/delx(i+1)
     columns(innz)=l-3
   else 
     innz=innz+1
     nnzval(innz)=1/delx(i+1)*2*eta(i,j)/(delz(j)+delz(j+1))        &  !u(i,j-1/2) =
                  -cpr*2/(delz(j)+delz(j+1))*2d0/3d0*etaS/delx(i+1)
     columns(innz)=l-4
   endif

   if(j.ne.1) then
     innz=innz+1
     nnzval(innz)=2/(delz(j)+delz(j+1))*2*etaS/delz(j)            & !v(i+1/2,j-1) 
                  -cpr*2/(delz(j)+delz(j+1))*2d0/3d0*etaS/delz(j)
     columns(innz)=l-3
   endif 

   innz=innz+1
   nnzval(innz)=-1/delx(i+1)*2*eta(i,j)/(delz(j)+delz(j+1))       & !u(i,j+1/2)   
                +cpr*2/(delz(j)+delz(j+1))*2d0/3d0*etaN/delx(i+1)
   columns(innz)=l-1

   innz=innz+1
   if(i.ne.nx) then
   nnzval(innz)=-2/(delz(j)+delz(j+1))*(2*etaS/delz(j)+2*etaN/delz(j+1))         & !v(i+1/2,j)
                -1/delx(i+1)*2*eta(i+1,j)/(delx(i+1)+delx(i+2))                  &
                -1/delx(i+1)*2*eta(i,j)/(delx(i)+delx(i+1))                      &
                +cpr*2/(delz(j)+delz(j+1))*2d0/3d0*(etaN/delz(j+1)+etaS/delz(j))
   else
   nnzval(innz)=-2/(delz(j)+delz(j+1))*(2*etaS/delz(j)+2*etaN/delz(j+1))         & !v(i+1/2,j)
                -1/delx(i+1)*2*eta(i,j)/(delx(i)+delx(i+1))                      &
                -1/delx(i+1)*bcR(j)*2*eta(i+1,j)*2/(delx(i+1)+delx(i+1))         &
                +cpr*2/(delz(j)+delz(j+1))*2d0/3d0*(etaN/delz(j+1)+etaS/delz(j))
   endif
   columns(innz)=l

   if(j.ne.nz) then
     innz=innz+1
     nnzval(innz)=2/(delz(j)+delz(j+1))*2*etaN/delz(j+1)            &  !v(i+1/2,j+1)
                  -cpr*2/(delz(j)+delz(j+1))*2d0/3d0*etaN/delz(j+1)
     columns(innz)=l+3
   endif
   
   if(i.ne.nx) then

     if(i.eq.1) then

       if(j.eq.1) then
         innz=innz+1
         nnzval(innz)=-1/delx(i+1)*2*eta(i+1,j)/(delz(j)+delz(j+1))       &  !u(i+1,j-1/2) * 
                      +cpr*2/(delz(j)+delz(j+1))*2d0/3d0*etaS/delx(i+1)
         columns(innz)=l+(3*nz-2)
       else
         innz=innz+1
         nnzval(innz)=-1/delx(i+1)*2*eta(i+1,j)/(delz(j)+delz(j+1))       &  !u(i+1,j-1/2) * 
                      +cpr*2/(delz(j)+delz(j+1))*2d0/3d0*etaS/delx(i+1)
         columns(innz)=l+(3*nz-3)
       endif
       
       innz=innz+1
       nnzval(innz)=2/(delz(j)+delz(j+1))-ala*Di*rho_r/2d0                   !p(i+1/2,j-1/2) +
       columns(innz)=l+(3*nz-1)
       
       innz=innz+1
       nnzval(innz)=1/delx(i+1)*2*eta(i+1,j)/(delz(j)+delz(j+1))       &     !u(i+1,j+1/2) -
                    -cpr*2/(delz(j)+delz(j+1))*2d0/3d0*etaN/delx(i+1)
       columns(innz)=l+(3*nz)
       
       innz=innz+1
       nnzval(innz)=1/delx(i+1)*2*eta(i+1,j)/(delx(i+1)+delx(i+2))           !v(i+3/2,j) #
       columns(innz)=l+(3*nz+1)
       
       innz=innz+1
       nnzval(innz)=-2/(delz(j)+delz(j+1))-ala*Di*rho_r/2d0                  !p(i+1/2,j+1/2) =
       columns(innz)=l+(3*nz+2)
       
     else

       if(j.eq.1) then
         innz=innz+1
         nnzval(innz)=-1/delx(i+1)*2*eta(i+1,j)/(delz(j)+delz(j+1))     &    !u(i+1,j-1/2) * 
                      +cpr*2/(delz(j)+delz(j+1))*2d0/3d0*etaS/delx(i+1)
         columns(innz)=l+(3*nz-1)
       else
         innz=innz+1
         nnzval(innz)=-1/delx(i+1)*2*eta(i+1,j)/(delz(j)+delz(j+1))     &    !u(i+1,j-1/2) * 
                      +cpr*2/(delz(j)+delz(j+1))*2d0/3d0*etaS/delx(i+1)
         columns(innz)=l+(3*nz-2)
       endif
       
       innz=innz+1
       nnzval(innz)=2/(delz(j)+delz(j+1))-ala*Di*rho_r/2d0                   !p(i+1/2,j-1/2) +
       columns(innz)=l+(3*nz)
       
       innz=innz+1
       nnzval(innz)=1/delx(i+1)*2*eta(i+1,j)/(delz(j)+delz(j+1))      &      !u(i+1,j+1/2) -
                    -cpr*2/(delz(j)+delz(j+1))*2d0/3d0*etaN/delx(i+1)
       columns(innz)=l+(3*nz+1)
       
       innz=innz+1
       nnzval(innz)=1/delx(i+1)*2*eta(i+1,j)/(delx(i+1)+delx(i+2))           !v(i+3/2,j) #
       columns(innz)=l+(3*nz+2)
       
       innz=innz+1
       nnzval(innz)=-2/(delz(j)+delz(j+1))-ala*Di*rho_r/2d0                  !p(i+1/2,j+1/2) =
       columns(innz)=l+(3*nz+3)
       
     endif
     
   else  
     innz=innz+1
     nnzval(innz)=2/(delz(j)+delz(j+1))-ala*Di*rho_r/2d0                     !p(i+1/2,j-1/2) +
     columns(innz)=l+(3*nz+1-2*j)
     
     innz=innz+1
     nnzval(innz)=-2/(delz(j)+delz(j+1))-ala*Di*rho_r/2d0                    !p(i+1/2,j+1/2) =
     columns(innz)=l+(3*nz+2-2*j)
   endif

   rhs(l)=loadz+kE*bcR(j)*1/delx(i+1)*(-2*eta(i+1,j)*2/(delx(i+1)+delx(i+1))*vR(j))

 end subroutine mom_v

!================================================
 subroutine cont_p(i,j,l,innz)
! Continuity equation for nodes p(i-1/2,j-1/2)
!================================================
   use share; use domain
   implicit none
   integer, intent(in) :: i,j
   integer, intent(inout) :: l,innz
   real (kind=8) :: z,rho_r,drho_r,T_r
!================================================

   if(cpr.eq.1) then
     z=zgrid(j)+delz(j+1)/2
     call refstate(z,rho_r,drho_r,T_r)
   else
     rho_r=1d0
     drho_r=1d0
     T_r=0d0
   endif

!  NB: The row built in this subroutine has zero element on the diagonal
!      that must be stored anyway (since symmetric PARDISO requires all diagonal elements)

   l=l+1
   rowind(l)=innz+1

!   if(rgx.eq.0.or.rgz.eq.0) then
   if(j.eq.0) then  
     
     if(i.ne.1) then

       if(i.eq.2) then
         innz=innz+1
         nnzval(innz)=-rho_r*1/delx(i)        !u(i-1,j-1/2) *
         columns(innz)=l-(3*nz+2)
       else
         innz=innz+1
         nnzval(innz)=-rho_r*1/delx(i)        !u(i-1,j-1/2) *
         columns(innz)=l-(3*nz+3)
       endif

     endif

     if(i.eq.1) then
       innz=innz+1
       nnzval(innz)=rho_r*1/delz(j+1)    &       !v(i-1/2,j) +
                    +cpr*drho_r/2  
       columns(innz)=l-(nz+1)
     endif

     if(i.eq.2) then
       innz=innz+1
       nnzval(innz)=rho_r*1/delz(j+1)    &       !v(i-1/2,j) +
                    +cpr*drho_r/2  
       columns(innz)=l-3*nz+1
     endif

     if(i.ne.1.and.i.ne.2) then
       innz=innz+1
       nnzval(innz)=rho_r*1/delz(j+1)    &      !v(i-1/2,j) +
                    +cpr*drho_r/2  
       columns(innz)=l-3*nz
     endif

     innz=innz+1
     nnzval(innz)=rho_r*1/delx(i)               !u(i,j-1/2)
     columns(innz)=l-1
   
   else 
     
     if(i.eq.1) then

       innz=innz+1
       nnzval(innz)=-rho_r*1/delz(j+1)   &      !v(i-1/2,j-1)
                    +cpr*drho_r/2  
       columns(innz)=l-(nz+2*j+2)

       innz=innz+1
       nnzval(innz)=rho_r*1/delz(j+1)    &      !v(i-1/2,j)
                    +cpr*drho_r/2  
       columns(innz)=l-(nz+2*j+1)

       innz=innz+1
       nnzval(innz)=rho_r*1/delx(i)             !u(i,j-1/2)
       columns(innz)=l-2

     elseif(i.eq.2) then 

       innz=innz+1
       nnzval(innz)=-rho_r*1/delx(i)              !u(i-1,j-1/2)
       columns(innz)=l-(3*nz+3)

       innz=innz+1
       nnzval(innz)=-rho_r*1/delz(j+1)    &       !v(i-1/2,j-1)
                    +cpr*drho_r/2  
       columns(innz)=l-(3*nz+2)

       if(j.ne.nz) then
         innz=innz+1
         nnzval(innz)=rho_r*1/delz(j+1)   &       !v(i-1/2,j)
                    +cpr*drho_r/2  
         columns(innz)=l-(3*nz-1)
       endif

       innz=innz+1
       nnzval(innz)=rho_r*1/delx(i)               !u(i,j-1/2)
       columns(innz)=l-2

     else

       innz=innz+1
       nnzval(innz)=-rho_r*1/delx(i)          !u(i-1,j-1/2)
       columns(innz)=l-(3*nz+4)

       innz=innz+1
       nnzval(innz)=-rho_r*1/delz(j+1)   &    !v(i-1/2,j-1)
                    +cpr*drho_r/2  
       columns(innz)=l-(3*nz+3)

       if(j.ne.nz) then
         innz=innz+1
         nnzval(innz)=rho_r*1/delz(j+1)    &    !v(i-1/2,j)
                      +cpr*drho_r/2  
         columns(innz)=l-3*nz
       endif

       innz=innz+1
       nnzval(innz)=rho_r*1/delx(i)                  !u(i,j-1/2)
       columns(innz)=l-2

     endif

   endif

!   endif

   innz=innz+1
   nnzval(innz)=0d0
   columns(innz)=l

 end subroutine cont_p


!=====================================================================
 subroutine refstate(z,rho_r,drho_r,T_r)
! Given z returns the value of the depth-dependent reference 
! density rho_r(z), its z-derivative drho_r(z) and the adiabatic
! temperature gradient T_r
!=====================================================================
   use share; use domain
   implicit none
   real (kind=8) :: z,rho_r,drho_r,T_r
!=====================================================================

! Constant expansivity
   if(alpha_law.eq.0) then
     rho_r = dexp(Di*(dimz-z))
     drho_r = -Di*rho_r
     T_r = Tsurf*(dexp(Di*(dimz-z)-1d0))

! Variable expansivity (but constant among different phases)
   elseif(alpha_law.eq.1) then
     rho_r = (alpha3*Di*(dimz-z)+1d0)**(1d0/alpha3) 
     drho_r = Di*(alpha3*Di*(dimz-z)+1d0)**(1d0/alpha3-1d0)
     T_r = Tsurf*(alpha3*Di*(dimz-z)+1d0)**(1d0/alpha3)

! Variable expansivity (with variable coefficients for different phases)
!   elseif(nphases.ne.0.and.alpha_law.eq.2) then

!     if(nphases.ne.0) then
!       do k=1,nphases
!         z0=phase_z(k)
!         Cl=phase_Cl(k)
!         Tph=phase_T(k)
!         ff(k)=fph(z)
!       end do
!     endif
!
!     if((dimz-z).le.z0_ol_sp) then
!       rho_r = (alpha3_ol*Di*(dimz-z)+1d0)**(1d0/alpha3_ol) 
!       drho_r = Di*(alpha3_ol*Di*(dimz-z)+1d0)**(1d0/alpha3_ol-1d0)
!       T_r = Tsurf*(alpha3_ol*Di*(dimz-z)+1d0)**(1d0/alpha3_ol)
!
!     elseif((dimz-z).gt.z0_ol_sp.and.(dimz-z).le.z0_sp_pv) then
!       rho_r = (alpha3_sp*Di*(dimz-z)+1d0)**(1d0/alpha3_sp) 
!       drho_r = Di*(alpha3_sp*Di*(dimz-z)+1d0)**(1d0/alpha3_sp-1d0)
!       T_r = Tsurf*(alpha3_sp*Di*(dimz-z)+1d0)**(1d0/alpha3_sp)
!
!     elseif((dimz-z).gt.z0_sp_pv) then
!       rho_r = (alpha3_pv*Di*(dimz-z)+1d0)**(1d0/alpha3_pv) 
!       drho_r = Di*(alpha3_pv*Di*(dimz-z)+1d0)**(1d0/alpha3_pv-1d0)
!       T_r = Tsurf*(alpha3_pv*Di*(dimz-z)+1d0)**(1d0/alpha3_pv)
!
!     endif

   endif


 end subroutine


!=====================================================================
 subroutine writesol
! Write the two components of the flow at non-staggered nodes 
! and pressure at staggered nodes
!=====================================================================
   use share; use domain
   implicit none
   integer :: i,j,ip
   real (kind=8) :: x,x1,x2,z,z1,z2,vx,vx1,vx2,vz,vz1,vz2,etam
!================================================================

   do i=0,nx+1
     x=xgrid(i) 
     if(i.eq.0) then
       x1=delx(1)/2d0
       x2=xgrid(1)+delx(2)/2d0
     elseif(i.eq.nx+1) then
       x1=xgrid(nx)-delx(nx)/2d0
       x2=xgrid(nx)+delx(nx+1)/2d0
     else
       x1=x-delx(i)/2d0
       x2=x+delx(i+1)/2d0
     endif
!     if(i.ne.0)    x1=x-delx(i)/2d0
!     if(i.ne.nx+1) x2=x+delx(i+1)/2d0

     do j=0,nz+1   
       z=zgrid(j) 
       if(j.eq.0) then
         z1=delz(1)/2d0
         z2=zgrid(1)+delz(2)/2d0
       elseif(j.eq.nz+1) then
         z1=zgrid(nz)-delz(nz)/2d0
         z2=zgrid(nz)+delz(nz+1)/2d0
       else
         z1=z-delz(j)/2d0
         z2=z+delz(j+1)/2d0
       endif
!       if(j.ne.0)    z1=z-delz(j)/2d0
!       if(j.ne.nz+1) z2=z+delz(j+1)/2d0
       
!-------------------------------------------------------------------
! Horizontal flow (averaged at Eulerian non staggered nodes)
! NB at the top (j=nz+1) and bottom (j=0) boundaries the 
!    velocity is extrapolated from the two closest solution points
!-------------------------------------------------------------------
!       if((i.eq.0.and.j.eq.0).or.(i.eq.0.and.j.eq.nz+1).or.  &
!          (i.eq.nx+1.and.j.eq.0).or.(i.eq.nx+1.and.j.eq.nz+1)) then
!         vx=0d0
       if(i.eq.0.or.i.eq.nx+1) then
         vx=0d0
       elseif(j.eq.0) then         
!         vx=upom((nz+1)*(i-1)+1)
         vx1=upom((nz+1)*(i-1)+1)
         vx2=upom((nz+1)*(i-1)+2)
         vx=vx2+(vx1-vx2)*(z-z2)/(z1-z2)
       elseif(j.eq.nz+1) then
!         vx=upom((nz+1)*i)
         vx1=upom((nz+1)*i-1)
         vx2=upom((nz+1)*i)
         vx=vx2+(vx1-vx2)*(z-z2)/(z1-z2)
       else
         vx1=upom((nz+1)*(i-1)+j) 
         vx2=upom((nz+1)*(i-1)+j+1) 
         vx=vx1+(vx2-vx1)*(z-z1)/(z2-z1)
       endif
       u(i,j)=vx       

!-------------------------------------------------------------------
! Vertical flow (averaged at non staggered nodes)
! NB at the left (i=0) and right (i=nx+1) boundaries the 
!    velocity is extrapolated from the two closest solution points
!-------------------------------------------------------------------
!       if((j.eq.0.and.i.eq.0).or.(j.eq.0.and.i.eq.nx+1).or.         &
!          (j.eq.nz+1.and.i.eq.0).or.(j.eq.nz+1.and.i.eq.nx+1)) then
!         vz=0d0
       if(j.eq.0.or.j.eq.nz+1) then
         vz=0d0
       elseif(i.eq.0) then
!         vz=vpom(j)
         vz1=vpom(j)
         vz2=vpom(j+nz)
         vz=vz2+(vz1-vz2)*(x-x2)/(x1-x2)
       elseif(i.eq.nx+1) then
!         vz=vpom(nz*nx+j)
         vz1=vpom(nz*nx+j-nz)
         vz2=vpom(nz*nx+j)
         vz=vz2+(vz1-vz2)*(x-x2)/(x1-x2)
       else
         vz1=vpom(nz*(i-1)+j)
         vz2=vpom(nz*(i-1)+j+nz)
         vz=vz1+(vz2-vz1)*(x-x1)/(x2-x1)
       endif
       v(i,j)=vz

     end do
   end do

!------------------------------------------------------------------------------------------------
! For non-linear models, correct the current iterate by applying a relaxation factor if necessary
! By default relax_fac=1 and no relaxation is applied
!------------------------------------------------------------------------------------------------
   if(nlv.eq.1) then
     u=relax_fac*u + (1d0-relax_fac)*u_previt
     v=relax_fac*v + (1d0-relax_fac)*v_previt
   endif

!---------------
! RMS velocity
!---------------
   do i=0,nx+1
     do j=0,nz+1
       urms(i,j)=dsqrt(u(i,j)**2 + v(i,j)**2)
     end do
   end do

!--------------------------------------------
! Various flow derivatives at (i,j) nodes
!--------------------------------------------
   do i=0,nx+1
     do j=0,nz+1
       
       if(i.eq.0) then
         dudx(i,j)=(u(1,j)-u(0,j))/delx(1)
       elseif(i.eq.nx+1) then
         dudx(i,j)=(u(nx+1,j)-u(nx,j))/delx(nx+1)
       else
         dudx(i,j)=(u(i+1,j)-u(i-1,j))/(delx(i)+delx(i+1))
       endif

       if(j.eq.0) then
         dvdz(i,j)=(v(i,1)-v(i,0))/delz(1)
       elseif(j.eq.nz+1) then
         dvdz(i,j)=(v(i,nz+1)-v(i,nz))/delz(nz+1)
       else
         dvdz(i,j)=(v(i,j+1)-v(i,j-1))/(delz(j)+delz(j+1))
       endif

       if(j.eq.0) then
         dudz(i,j)=(u(i,1)-u(i,0))/delz(1)
       elseif(j.eq.nz+1) then
         dudz(i,j)=(u(i,nz+1)-u(i,nz))/delz(nz+1)
       else
         dudz(i,j)=(u(i,j+1)-u(i,j-1))/(delz(j)+delz(j+1))
       endif

       if(i.eq.0) then
         dvdx(i,j)=(v(1,j)-v(0,j))/delx(1)
       elseif(i.eq.nx+1) then
         dvdx(i,j)=(v(nx+1,j)-v(nx,j))/delx(nx+1)
       else
         dvdx(i,j)=(v(i+1,j)-v(i-1,j))/(delx(i)+delx(i+1))
       endif

!    Viscous dissipation Phi = (sig_ij*eps_ij)/2  at (i,j) nodes
       Phi(i,j)=2*eta(i,j)*(dudx(i,j)**2+dvdz(i,j)**2)          &
                +eta(i,j)*(dudz(i,j)+dvdx(i,j))**2              &
                -cpr*2d0/3d0*eta(i,j)*(dudx(i,j)+dvdz(i,j))**2

!    Second invariant of the strain-rate		
!      eII=sqrt(e_ij e_ij)
       epsII(i,j)=dsqrt((dudx(i,j)-cpr*1d0/3d0*(dudx(i,j)+dvdz(i,j)))**2+0.5d0*(dudz(i,j)+dvdx(i,j))**2 &
                        +(dvdz(i,j)-cpr*1d0/3d0*(dudx(i,j)+dvdz(i,j)))**2)

     end do
   end do

!------------------------------------------------------------------
! Pressure and vertical deviatoric stress at staggered nodes
! etam: viscosity at p-nodes computed from neighboring 4 grid nodes
! vz1: vertical flow below the p-node
! vz2: vertical flow above the p-node
!-------------------------------------------------------------------
   ip=0
   do i=1,nx+1
     do j=1,nz+1
!  pressure node in the upper left corner is excluded from the solution
!  because of the condition of pressure uniqueness                
       if(i.ne.1.or.j.ne.nz+1) then 
         ip=ip+1                   
         p(i,j)=ppom(ip) 

!    Average viscosity at staggered nodes
         etam=(eta(i-1,j-1)+eta(i,j-1)+eta(i,j)+eta(i-1,j))/4
         
         if(j.eq.1) then 
           vz1=0d0
         else
           vz1=vpom(nz*(i-1)+j-1)
         endif
         if(j.eq.nz+1) then
           vz2=0d0
         else
           vz2=vpom(nz*(i-1)+j)
         endif
         
       endif

     end do
   end do

 end subroutine writesol


!#########################
 end module flow
!#########################

